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Abstract 

We define a set of operators that localise a radial image in radial space and radial frequency 
simultaneously. We find the eigenfunctions of this operator and thus define a non-separable orthogonal 
set of radial wavelet functions that may be considered optimally concentrated over a region of radial 
space and radial scale space, defined via a doublet of parameters. We give analytic forms to their energy 
concentration over this region. We show how the radial function localisation operator can be generalised 
to an operator, localising any L^(R^) function. We show that the latter operator, with an appropriate 
choice of localisation region, approximately has the same eigenfunctions as the radial operator. 
Based on the radial wavelets we define a set of quaternionic valued wavelet functions that can extract local 
orientation for discontinuous signals and both orientation and phase structure for oscillatory signals. The 
full set of quaternionic wavelet functions are component wise orthogonal; hence their statistical properties 
are tractable, and we give forms for the variability of the estimates of the local phase and orientation, 
as well as the local energy of the image. By averaging estimates across wavelets, a substantial reduction 
in the variance is achieved. 
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I. Introduction 



Localised analyses in one dimension have proven to be remarkably successful - notably so wavelet 
analyses. The latter is based on the idea that observed signals varying over an increasing argument, time 
say, exhibit disparate and highly localised behaviour associated with variations at a particular scale and 
at particular time points. Analysis is based on the wavelet transform, given for signal g{t) using wavelet 

1^(1) via 



where a is referred to as the scale, b the translation and * denotes conjugation. This allows for the 
recognition of patterns specific to time localisations b and length scales associated with scale o, if the 
function -(/'(•) is chosen such that the support of is essentially limited to a region near the origin, 
and the support of the Fourier transform of V'(")i is essentially limited to a region near some non-zero 
reference frequency /max- A function cannot be perfectly compact in time and frequency simultaneously, 
and so other criteria have been specified to determine the localisation properties of Of particular 
note is the idea of a localisation operator, generalizing the truncation in time or frequency operators 
[1] to simultaneously localising in time and frequency/scale [2], [3]. The eigenfunctions/eigenvectors of 
such operators are optimally localised with respect to the operator and in one dimension the problem of 
defining appropriate operators and calculating their eigenfunctions has been considered in detail [2], [3], 



The choice of extension of decomposition to two dimensional analysis is not trivial since variation in 
the spatial variable is associated with a direction, as well as a scale. This direction cannot be assumed 
to be aligned with the observational coordinate axes, and thus analysis using a simple tensor product of 
one-dimensional wavelets is, in general, not suitable. In two dimensions localisation is made to spatial 



w^{a,b; g) 



(1) 




[4]. 



point b = [bi, 62]^ , in scale to a and in orientation to angle 9 G [0,27r) cf [5]. The two dimensional 
continuous wavelet decomposition of image g{x) using wavelet tlj{x) is constructed via 



■w^{a,9,b; 



9) 



{i'a,e,b,9) 



(2) 




where (Da) represents a dilation, (Tf,) a translation and (TZg) a rotation of tpix) giving 



i'a,e,b (x) 



(3) 




(4) 
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with rg given as the rotation matrix. The decomposition will, with an appropriate choice of wavelet 
function, uncover/disentangle behaviour across specific spatial points, scales and orientation. Local two 
dimensional patterns in general may be intrinsically one dimensional, i.e. after a suitable rotation all 
variation is along a single axis, or intrinsically two dimensional, i.e. there is variation in several directions, 
operating at the same scale emanating from one spatial point. Following remarks by [6] we focus on 
wavelet analysis of discontinuities and oscillatory structure. Note that edges, or spatial discontinuities, 
have an orientation if they locally correspond to (one dimensional) curved discontinuities whilst (two 
dimensional) point discontinuities have no associated orientation. Oscillations may structurally take the 
form of one dimensional objects such as repeated lines with an even spacing that, if rotated to the 
appropriate axes, can be locally described as constant in one variable and as a sinusoid in the other. Two 
dimensional oscillations, circularly emanating from a single point, when considered locally at a distance 
from their source may be described approximately as one dimensional oscillations. 
Thus the structure of an image is highly orientation dependent, and analysis methods should disentangle 
both locally one dimensional and two dimensional structures operating at many different orientations. 
A well known feature of wavelet analysis, [7] is that genuinely two dimensional structure, i.e. point 
discontinuities, are well represented in a wavelet decomposition, and so in this paper we focus instead 
on the treatment of locally one dimensional structures by adjusting the wavelet transform suitably. We 
shall discuss two existing strategies for considering oriented scale -based decompositions, and construct 
a new method corresponding to a synthesis of the two methods discussed. This method can extract the 
local orientation of the image explicitly, in a multi-scale framework. 

Existing continuous wavelet methods that deal with the orientation of the image explicitly are based on 
directionally selective filters, or directional wavelets. Antoine & Murenzi [5] define complex directional 
wavelets with a preferred orientation in the frequency domain, as their frequency support is limited to a 
pre-defined cone, parameterised via the opening and closing angles of the cone [5, p. 324-5]. Defining 
highly directional wavelets will necessitate an elongation of the wavelets in the spatial frequency domain, 
and this affects their spatial and spatial frequency resolution capacity - along a specific orientation in 
the frequency domain, the wavelets localise badly in frequency, and for this reason we do not use 
directional wavelets. Directional wavelets will localise in scale and direction, whilst spatially isotropic 
wavelets only separate features at different scales. We shall use an isotropic wavelet decomposition to 
separate out disparate components occurring at either different scales and the same spatial locations, 
or at the same scales at different spatial locations. To facilitate this separation of structure, wavelets 
that are optimally concentrated in radial space and radial frequency are required, as in two dimensions 
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the notion of spatial distance is naturally associated with the Cartesian metric. We define a family of 
radial two-dimensional localisation operators and find the radial eigenfunctions of any given operator 
in this family, denoted the isotropic Morse wavelets. Any operator in this family is characterised via 
two parameters that determine the spatial/spatial frequency structure of the isotropic eigenfunctions. Any 
choice of the parameters fixes a particular operator that in turn possesses a family of eigenfunctions. These 
functions are orthogonal and indexed via n. The eigenvalues explicitly give the radial concentration of 
the eigenfunctions. These eigenfunctions are related (but not equivalent) to the eigenfunctions of the one 
dimensional Morse localisation operator [3]. For every fixed value of n, and given radial eigenfunction, we 
define an additional pair of functions, whose joint norm may be considered to have the same localisation 
in space and spatial scale as the original radial function, but when combined with the original radial 
function will disentangle the local orientation of the image analysed. These extra pairs of functions are 
constructed explicitly to consider local orientation and phase. 

A method for considering local phase structure is to extend the notion of instantaneous frequency [8]; 
[9], [10], [11] give extensions to instantaneous frequency and local phase structure in the spatial domain, 
and for each spatial point retrieve a local phase/variational structure. These extensions correspond to the 
calculation of several additional images, or quadrature components, at each spatial point, where each 
additional set of components is considered to have the same local spatial energy and variational structure 
as the original image. The full set of components is used to calculate the local orientation and variational 
structure. The additional components are ill-defined when constructed from multi-component images, 
as then a single component with a spatially varying phase function is not an appropriate model for 
the original image. It thus becomes necessary to combine the calculation of a local phase with scale- 
localised methods such as the wavelet transform. For each fixed value n G N, for each Morse wavelet, 
we define two extra real functions to complement the isotropic Morse wavelets. The triplet of real valued 
functions form a monogenic signal [9]. Two of the triplet of functions should be thought of as a single 
vector valued object, where their vector structure characterises the orientation of the local variations. The 
monogenic wavelets and wavelet transform are best represented using quaternion [12], rather than, real 
or complex numbers. Each triplet is therefore considered as a positive real valued amplitude, a pure unit 
quaternion specifying a direction, and a phase. The real amplitude characterises a local energy, the pure 
unit quaternion an orientation, and the phase a local variational structure [9]. The quaternion algebra 
allows for easy parameterisation of phase and orientation structure. 

In one dimension analytic wavelets, i.e. Cauchy or Morlet wavelets [13, p. 28], are used to identify local 
oscillatory structure, of a real image g{-). The monogenic Morse wavelets are the natural two-dimensional 
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extension of the analytic Morse wavelets [3], and they define a local phase and orientation structure at each 
spatial/spatial scale point for a two dimensional oscillatory image, yielding a natural, and more elegant, 
structure for wavelet ridge analysis [14]. Oscillatory images are the complements of discontinuous images, 
where the two kinds of images appear as either oscillatory or discontinuous depending on if they are 
analysed in the Fourier or spatial domain. We consider line, or curved, discontinuities that may locally 
be thought of as one dimensional structures, and discuss the retrieval of their local features. Analyses 
using the separable discrete wavelet transform in addition with local phase structure characterisations, 
have been previously considered [15], [16]. However, in contrast to their procedure, we define multiple 
orthogonal continuous wavelets, based on a different two dimensional extension to the analytic signal. 
The Morse wavelets are additionally optimally concentrated with respect to a radial position/scale region 
V. Other quaternion valued decompositions includes the work of [17], however this decomposition is 
only suitable for deterministic images. 

We briefly discuss the discrete implementation of the two dimensional monogenic Morse wavelet trans- 
form, and in more depth the statistical properties of the transform. Unavoidably, most observed signals 
are contaminated by noise, and so robust methods that can deal with noise must be designed. As the 
operator problem yields solutions of multiple orthogonal monogenic wavelets, we may use the notion 
of averaging uncorrected estimates [1] to retrieve estimates of the space-scale energy as well as other 
quantities of the image with reduced variance. Multiple orthogonal filters have been considered in several 
dimensions for stationary processes [18], and non-stationary processes [19] using the windowed Fourier 
transform and tensor product windows, but our wavelets are in contrast to these methods orthogonal, 
non-separable monogenic wavelet functions. Finally, the methods are illustrated on typical examples, 
showing the power of the multiple monogenic Morse wavelets. 

II. Notation 

We denote the ID Fourier transform of g{-), G{f) = (exp(2jf7r/t), (7) , and the two dimensional Fourier 
transform of g{x) as G{f) = {exp{2jTrfx),g). We denote an arbitrary quaternion via q = qi + (?2* + 
qsj + 94^) where G R, i = 1, . . . 4 and note that P = = = ijk = —1, whilst ij = —ji = k, 
ik = —ki = —j, and finally jk = —kj = i. The algebra is non-Abelian. We additionally define the two 
dimensional Fourier transform in terms of any unit quaternion as Ge^{f) = (exp(2eg7r/a;), (7), so that 
the regular Fourier transform corresponds to Gj{f) = G{f). For more notes on quaternion algebra see 
[12]. We retain here only the briefest possible usage of the quaternion algebra, necessary for clarity of 
exposition, and stress that all implementation is discussed in terms of real vector quantities. The rotation 
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Fig. 1. The isotropic Morse wavelets in the spatial domain. {l,m) — (8,3) and n = (far left), n = 1 (second from left). 
The modulus of the isotropic Morse wavelets in the spatial frequency domain. {l,m) — (8, 3) and n = (second from right), 
n — 1 (far right). 



operation is implemented using matrix 

_ /cos(0) -sin(e)\ 
\sm{e) cos{e) ) 

We shall extensively use polar coordinates, and define: x = [x cos{x),x sin(x)] , f = [f cos{(j)), f sm{(j))] , 
and b = [6cos((/)b), 6sin((/)fe)] . 

III. Two DIMENSIONAL WAVELET ANALYSIS 

One dimensional local analysis coiTcsponds to decomposing a function in terms of a set of functions that 
contain behaviour local to a set scale a, and time point b. The two dimensional wavelet transform is defined 
using four parameters: ^ = [a, 9, b] , where a and b play roughly the same role as the corresponding 
one dimensional quantities, and 9 corresponds to local orientation localisation. For the decomposition to 
be meaningful, the translated and dilated wavelets are chosen in one dimension to be mainly supported 
near time point b and frequency point fa = fma.x/a, and the obvious extension to two dimensions 
would be to find functions that are mainly supported at spatial point b and spatial frequency point 
fa = [cos(0), sin(0)] . To measure the localisation of a given wavelet function, commonly its spread 
in time and frequency or scale is calculated, and both quantities are desired to be low - their product 
is bounded below, thus restricting possible joint localisation. A function's spread in time and frequency 
should be considered simultaneously in the two domains [3], rather than combining two separate marginal 
properties: for this purpose localisation operators were defined, denoted Vx>- The localisation of a function 
g{-) to region V is measured by the ratio of energy figiV) = {V-Dg{-),V-Dg{-)) /{g{-), g{-)), and the 
eigenfunctions of Vx> achieve optimal ratios [3]. Naturally extending analysis to two dimensions requires 
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the appropriate definition of a two-dimensional localisation operator. We consider localisation in scale and 
spatial location, which requires localisation in radial spatial frequencies of spatially radial functions, as we 
wish to associate a notion of distance to the Cartesian two dimensional metric. This produces functions 
optimally concentrated in a radial space and scale a, but that have no orientation. We then construct a set 
of functions capable of extracting local orientation information based on the monogenic signal [9]. [20] 
have shown that with an appropriate definition of scale, orientation and spatial spread, the orientation and 
scale may be considered separately when finding optimally concentrated functions. In a slightly different 
setting they demonstrate that radial functions have optimal scale versus position localisation properties, 
for appropriately chosen families of functions. This, in combination with a Cartesian metric in the plane, 
motivates the study of radial localisation. 

A. Radial Localisation Operator 

The construction of a coherent radial state, required for the construction of a radial projection operator, 
is not straightforward. We commence with radial function v{x) = Vr{x) and in theory wish to construct 
a family of radial functions that have been shifted in scale and position. Obviously this is impossible, as 
once we shift in position, we no longer have a radial function, but we may relax our requirements as the 
family only needs to act in a similar fashion to a direction averaged spatial shift, based on an appropriate 
domain V chosen. To construct a generic two-dimensional non-radial coherent state we would use the full 
set of parameters ^ = [a, 9, h] , but intend to use the sub-set = [a, b] . Working with radial functions 
the act of rotation, represented by 6 will have little importance, and is not included. The dilation by a will 
act correspondingly in the world of radial operations, to one dimensional analysis, and a substitute for 
a shift in position by h must be defined. We reconstruct the function using building blocks of {v^,g)v^ 
that are averaged over suitably defined domains. The localisation domain will be defined in terms of b 
and range over < < 27r. This implies that any shift in b will be averaged over the full range of 
(j)b and so rather than multiplying V{f) by e~^^^f'^ we could multiply this by its angular average of 
Jo(27r/6), that after implementing the operator projection, yields the same results as the former strategy. 
Unfortunately, this choice of coherent state leads to a mathematically intractable operator, and so we 
define a radial coherent state for a > and 6 > that has similar properties to the suggested state via 




cos {uj'^b — vr/4) 



(5) 




based on the Morse coherent state [3] of 



2r/2+l/2^ 



w^e"^' if cj > 0. 
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Note that y^'T(tj) (X Lo^, (3 > for Lij << 1, and for future reference define u>^ : V^''^{lo) < 
e \/lu < lu^. In equation (5) the dilation of a is implemented as in [3] and requires no further discussion. 
The factor {a^/"'u))~^/'^^"'/'^ is added to ensure the correct normalisation of the two-dimensional radial 
function, as is the replacement of a^^"' for the term in a multiplying the dilated function of oj. 

Denote the translation-like operator, and the originally posited direction averaged operator respectively 
as, tI^^\lo, b) = ^nd Tlf^{uj, h) = Jo{uj'^b). The translation operator's decay for large values 

(2) 

of Lo'^b is the same as as is its zero crossing structure, and apart from small arguments the 



two functions are performing a similar action. The functional behaviour for small values of uj is different 

1 

V2uj-'b 



- Tf V: b) = l + (^2^) whilst tI'^\uj, b) = + O [uj^''"^) , where the latter is unbounded near 



a; —> 0. If we denote the coherent state as 

= V2ai/^(ai/^a;)-V2+^/V^'^(ai/^w)r,(")(cc.,6), n= 1,2, 
then for small values of with r = (2/3 + l)/7, we find 

so that both choices give contributions of negligible magnitude for — > as long as we assume 
P > 1/2, which is combined with the previous constraint of 7 > 1, with /3 > (7 — l)/2 [3, p. 2663]. We 
choose T^^^^ as the radial spatial shift quantity. This has an approximate interpretation of 6 as a orientation 
averaged spatial shift. Note that the generalised Morse wavelets also are based on the approximate, rather 
than exact, notion of a warped location shift [3, p. 2663]. 

We define the operator for a radial function g{x) in terms of radial inner product {gi,g2)R = 
^{Gi,G2)r = ! G\{f)G2{f)f df as and (in terms ofV+{C) = {(a, 6) : a'^ + b'^ + l< 2aC,b> O}), 

da 



db 



Go / / K2{uJi,uJ2;a,b) ^ db du2, (6) 

Juj2=0 J Ja^+b^+l<2aC « 



where 



K2{uJi,u;2;a,b) = 2a^/Ty^'^(ai/^a;i)vW^'t'^'^(«'/^^2)G(^)cos((cj7 - u^)b). (7) 

zvr 

The kernel K2(-, ■;■,■) can be inverted in its first argument to retrieve the spatial domain operator. Note 
that by definition cui > 0. We constrain the norm of the operator - exactly reconstructing the entire radial 
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function if we let the region of integration across be large enough, thus calibrating the operator to 
make the eigenvalues meaningful. 



. poo poo poo 

Juj2=0 Ja=0 Jb=0 



OO POO POO 



J uj2=0 J a=0 Jb=—oo 



[^) cos{{ujJ - uj'^)b) ^ 



= Co2(27r)2G(^)^, (8) 
2-ir r — I 

which gives Co = yielding a 'resolution of identity' [3] for radial functions: this will not hold 

for any g{-) G L^(R^), as the operator is only defined for radial functions. The eigenfunctions of the 
operator defined in (6) can be found by solving the equation 

^^^WG}(|1) = AG(^). (9) 

The one-dimensional Morse wavelet projection operator can be considered to be in the case of tpit) real, 
assuming that toi > 0, (similar expressions are derived for cji < 0, but as we shall use this to obtain 
solutions to equation (9) we only need to consider cui > 0, and the sin(-) term vanishes due to the 
symmetry of the projection region, and where V is as in [3]): 

= Cl^ r [ [ a^/''V'^'^{a^/"'ui)V'^'^*{a^/''uj2) (10) 

^2tt^ ^ ^ 1 ^2 27r 

tlD roo 



where 



[ II Hii{uJi, u;2, a, b) db du;2, 

Juj2=0J Ja^+b^+l<2aC, aeR+,beR+ « 



Ki(wi,u;2;a,6) =2ai/^y'3'^(ai/^wi)y^'T(a^/^u;2)^(— )cos((u;7 (11) 



Note that the above kernel •; •, •) is similar to the kernel K2(-, •; •, •) of (7), the only difference being 
that K2{-, ■;■,■) has the extra term, ^u)2/^i, and = The Morse wavelets [3] are the solution 
to the equation 

T'h''mi^) = X'^n^). (12) 
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Consider equation (9); multiply both sides by ycJi7(27r), set G{uJ)^/uJ = "^{to), and note that the equation 
to be solved has now exactly the form of equation (12): 

where we have used that 27rCo/C^^ = 1. Thus the solutions of (9) are given by 

^'r!^l,M) = 4f'^iLif), f>0, (14) 



where ^'^^l (•) are the even Morse wavelets as defined in one dimension [3], / = /? — |, m = 7, and 



n ^ N enumerates the eigenvectors. The eigenvalues correspond to, 



C-l 



r(n + i)r(r-i) 70 

and this yields the radial-spatial, radial-scale concentration of ('D(C)) = A^r(C')- It may seem 
surprising that, in two-dimensions, the same eigenvalues, and thus concentration values, are found as in 
the one-dimensional case, note, however, that this is only derived for radial images, that are constrained 
to the same behaviour in both spatial directions - hence in essence we are really only making a one 
dimensional compromise. The hypervolume of V is directly related to C, and note that we may formulate 
the notion of bias of estimation of local properties of the signal, or leakage, in terms of the eigenvalues 
as 1 - Xl^C) [3]. 

B. Non-Radial Localisation Operators 

The operator outlined above was constructed in the radial frequency domain for explicitly radial images. 
Let us explicitly consider how a operator is constructed in two dimensions. Define 

V2 = {$:a'^ + b'^ + l<aC, 0<e <2tt, b G a > O} . (16) 

Let us again consider a coherent state v'^{x) that is the building block of the projection operator, but let 
us now make this local to ^. Define the Fourier transform of this coherent state at ^ with 2nf = u as 

V+{u) = al/(27)+l/2^-l/2+7/2y/3,7(al/7^)e-jb{r--.a,V--^ 

= ^ilj^!i^V^a{r+l)/2^-l/2+7/2^/3g-aa;^-j5{r_«a;)a;^-i ^-^j^ 

r(r) 

The normalization of V~^, for completeness, can be calculated by: 

{v^.,v'^) = a^^^~^^ ^^1)2 Jq^ uj~^~^'^V^'^'^{a^^^uj)u!dujd(f) = |. Therefore, the coherent states are of 
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norm unity, if they are multiplied by y^. We define the localisation operator for any function g{x) E 
L2(m2) as 

Vv, {g} (x) = Coj^ v^{x){vl,g) ^SbdO. (18) 

Equation (18) gives an expression for the localisation of an arbitrary function g{-) over region X>2- We can 
by calculating ^g{V2) find the localisation of gQ G L^(M^), to ■ We shall take Co such that as ID2I — > 
oo> 'Pv2 {9} (^) ~^ 9{x) for ^11 radially symmetric functions. Note that if g{-) is a radially symmetric 
function, then {v+,g) = ^ aV(27)+i/2^-i/2+7/2y/3,7(ai/7cj)G(i^) Jo(cjT&)u; duj, where Jo(-) is 
the zeroth Bessel function [21]. Integration over (/>, for a radial g{-), thus removes the angular dependence 
on 6 and (pb, see [22]. In the frequency domain 

PvAG}M = Co j^v^{uji){v^,c)^d'hde 



00 



= 27ra / / / al/(27)+l/2^-V2+7/2^A,(^l/,^^)^l/2+7/2^1/(2,)+l/2 

xV^^^{a^''<U2)G(^)jQ{bujl)Mbujl) ^bdbdL02 (19) 

= Co [ [ [ K{uJi,uJ2;a,b) db du>2- (20) 

Note that V^''^{uj) ^ if uj < uj^ and thus we consider the term Jo(w^6) for ujjb » 0, as the point 
6 = has zero measure in the plane. For \z\ » the asymptotic approximation to the zeroth Bessel 
function is 

Jo{z) -^cos(z - 



cf [21, p. 364, 8.2.1]. We require |argz| < ir however, this is not an issue as 6 > and toi > 0. Thus 
for fixed non-zero b for values such that the integrand is non-zero 

K{ui,u;2;a,b) = 27raV7,,^-V2+7/2^A7(„l/7^^)^l/2+7/2^/3,7(„l/7^2)G^(^) 
UbLoDUbLoDb 

^ 27rai/^u;-^/'+^/V^'^(aV7^^)^V2+7/2^A,(^i/,(^^)^(^^ 

27r 



^ ■cos(a;^6- ^)J^y-cos(a;76- 



= 2ai/^F^'T(ai/^a;i) vW^^^'^(o^^^W2)G(^) cos((a;7 - a;^)5) 

+2ai/^y^'^(aVTc^i) vW^^'''^(a'^^(^2))G(?) cos((w7 + w^)6 - J). (21) 

zvr z 

The integration is over uj2 > and also ui > 0, hence the first term of (21) dominates over the second 
term. The integrand of k{loi,uj2', o,,b) can be replaced by ^2(^1, o, 6) defined in (7) of the previous 
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section. Thus we approximate the operator acting on G(-) in the frequency domain via 

~ f'°° f f da 
VvAG}{u^i) - Co // K2{uJi,uj2;a,b)^dbdL02 = Vv{G}{u:i), (22) 

Juj2=0 J Ja2+62+i<2aC 

and Co can be found from Co- Defining "Ppj allows for the consideration of the localisation of an arbitrary 
L^(M^) function. The radial eigenfunctions of 7^X12 are approximately those of Vx), where the derivation 
of the approximation shows the reasoning behind the definition of Vx>- Finally Vx>2 can be generalised 
to an arbitrary localisation by removing the constraint of radial symmetry in D2 and on g{-). 

C. Isotropic Wavelet Definition 

The even multiple Morse wavelets are defined in one dimension in the Fourier domain for fixed n = 
0,1,2..., /3>1, 7 > (/3-l)/2 and denoted by ^n-p,'r (■^^ of A„;^,^ = y^njlTin + l)/r(n + r) 

[3]. We define the isotropic two dimensional wavelets as the eigenfunctions of Vx>, for fixed n,l,m in 
terms of / = ||/|| via 

<;L (/) = ^(2vr/)'e-(2-/)'"L<- (2(2vr/)-) , (23) 
with c'l,^ = (2/ + 2)/?7i — 1, to attain the correct normalisation over 

L2(]R2) and where L^^ (•) are 

generalized Laguerre polynomials. The spatial domain wavelets with x = \\x\\ are found via the inverse 
Fourier transform for radial images, V^l^^ (x) = 27r (/) Jo (27r/a;) fdf. We plot the wavelets for 

/ = 8 and m = 3 with n = 0, 1 in the spatial domain, see Figure 1. Their radially symmetric oscillatory 
structure is very clear. They are optimally concentrated in a radial structure centred at the origin. We plot 
the modulus of the same function in the spatial frequency domain in the same plot. These functions are 
band-pass filters that are non-zero for a range of frequencies centred at the same distance from the origin 
in the frequency domain. The trough in the n = 1 follows as the first two wavelets are orthogonal, and 
we see that the sum of the moduli will be large in the same ring-shape structure. To further characterise 
the Morse wavelets define the radial frequency that maximises the isotropic Morse wavelets as 



/mix = argy>omax 



n = 0,...,N -I. (24) 



The magnitude square of the Fourier transform of analysis wavelet tp^^l^i-) will have a maximum at 
frequency fmlx/a and is unaffected by both b and the rotation. 

IV. Monogenic Images 

In one dimension the analytic signal is used to unambiguously define the phase and amplitude of a 
given real signal and using an analytic analysis wavelet will allow for the definition of a local magnitude 
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and phase at each time and scale point - a necessity for the analysis of multi-component images. The 
analytic signal is constructed in one dimension by removing all negative frequencies in the signal, and 
then inverting the Fourier transform - any real signal u{t) is complemented by its Hilbert transform 
v{t) = H{u}{t), and the analytic signal corresponds to u~^{t) = u{t) + iv{t). If an oscillation is 
persistent at a particular range of times, then it will be heavily weighted in the Fourier domain, and the 
analytic signal will approximately take the form of a complex exponential. 

The correct extension of the analytic signal to two dimensions has been the subject of much debate - 
of particular note are perhaps the single orthant image of [11], the hypercomplex signal of [10] and the 
monogenic signal of [9]. Following Felsberg and Sommer [9] we define the Riesz transform of an image 

u{x) as 

nu{x) = iniu{x) + j7^2n(a;) = iv^^\x) +jv^^\x). (25) 

The Fourier transforms of these two objects are 

^{7^ln(a;)} = -j cos (</.) [/(/), J-{7^2n(a;)} = -j sin (0) [/(/). 

Define the monogenic image [9] as u^'^^x) = u{x) +TZu{x). This is a quaternion valued object, and 
relations between the components of the quaternion are interpretable in terms of orientation and phase, 
for oscillatory images, as will be demonstrated. 

Consider an oscillation in two dimensions, corresponding to repeating maxima spaced l//o apait in 
orientation n = [cos (77) sin(ry)]^ . This corresponds to the simplest oscillatory image 

gi{x) = ai cos {lirf^x ■ n + Og) , (26) 

where ai, /o and 6s are constant scalars, whilst n is a constant unit length vector. Note that in variable 
y = X - n, gi{x) solves the one dimensional harmonic oscillator equation of ^ + (27r/o)^ 9i{y) = 0. 
To characterise a given image, we wish to determine ai, /o, 6s and n, from the image. 
We can calculate the monogenic extension of (26), and obtain, with e^i = icos(??) + jsin(7?), [23] 

g[-'\x) = g{x)+ig^'Hx)+jg^'Hx) 

= ai [cos {2Tr f qx ■ n + 6s) + erisin{2Tr /qx ■ n + 6s)] ■ (27) 

In quatemionic polar coordinates, 

g[-^\^) = aie^''''^^^"y+^^\ (28) 
and we may determine oi = \/ g'^{x) + (7(^)^(3;) + g^'^^'^{x), rj = tan"-*^ I 9!_H2Ll 



foy + 0s = ^tan^i ( 3(cc) ) • We restrict -| < r/ < §, and -i < /oy + 
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6s < ^. In terms of y we hence have a frequency domain description that perfectly mirrors the one 
dimensional theory, and retrieve the properties of the image from its monogenic extension. Naturally, in 
most real applications perfectly oscillatory images are not encountered, and more general models must 
be considered. 

V. AM/FM/OM Images 

As in [14] and [24], we shall consider images that locally may be approximated as a sum of sinusoidal 
components. This model is of some importance in machine vision, and can be applied to granular flow 
and general oriented patterns [25]. Assume that 

L 

c{x) = '^ci{x), ci{x) = ai{x)cos{2TT4>i{x)), (29) 
1=1 

where (t)i{x) = ipi{ni{x) • a;), with the added constraint that the unit vector ni{x) is varying slowly, in 
comparison with x, across the spatial period. We introduce this extra notation so that we may characterise 
images that can be considered as approximately sinusoidal in variable yi{x) = ni{x) ■ x. We refer to 
ni{x) as the orientation modulation (OM), whilst is the phase modulation of component /, and 
plays the same role as the phase/frequency modulation (FM) of a one dimensional signal. Then we find 
that for x = Xo + Sx, and ip'i{yi) = d/dyi [^i{yi)\ , that 

(t)i{x) = (piixo) + ipi{ni{xo) ■ Xo)ni{xo) ■ {x - Xo) + O (\Sx\'^) . (30) 



In the following we will use the shorthand ^^(a;) = ip'i{ni{x) ■ x). We additionally assume ai{x) varies 
slowly in comparison to the cosine term, and this corresponds to the amplitude modulation (AM). The 
monogenic version of the Ith component is 

cl^(^x) = a^(a;g)e^'^^'(^°K'^'(^°)+'^i('^°)'^'('"°)'(^~'"°)+*^(l^~'^°l^)) (31) 

= aiixo) [cos {2tt(/)i{xo) + 2TT(/)[{xo)ni{xo) • {x - Xo)) 

+ei{xo) sin {2n(l)i{xo) + 2'K<p'i{xo)ni{xo) • {x - Xo))] + o{\x - Xo\). 

c/ equation (27). Assuming the orientation ei{xo) to be stable across values of x for which ai{x) is non- 
zero, we shall perform the Fourier transform of c^~^\x) in terms of the unit quaternion ei{xo) instead 
of j. Note that de Moivre's theorem is still valid for any unit quaternion, and so this directional Fourier 
transform can be interpreted just like the regular Fourier transform in terms of oscillatory components. 
The directional Fourier transform is then given by 

C+^iif) = J J Mxo) e2^'='(^°)(<^'(^°)+'^;(^°)"'(^°)(^-^°))e-2^^'(^°)-'^^ d^x. (32) 
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h ^2 '? '2 

Fig. 2. The xi Riesz transform Morse wavelets in the spatial domain for (/, m) — (8, 3) and n — (far left), n = 1 (second 
left). The modulus of the xi Riesz transform Morse wavelets in the spatial frequency domain for (l,m) = (8,3) and n = 0, 
(second from right) n — 1 (far right). 

We apply the stationary phase approximation to this integral [14], as the unit quaternion will be acting 
like j itself. Thus, under the assumption that Q,{x,f) = 4>i{xo) + (l)'i{xo)ni{xo) ■ {x — Xq) — fx 
has the unique stationary point Xg, and there is quadratic behaviour around this point, i.e. ^}{x,f) = 
^ [{x — xq)H{x — xq)] , where H is the Hessian matrix of i}{x, /), we find at f{xo) = (p'i{xo)ni{xo), 
the integral provides the only non- null contribution of C+(/) W 2^^^(^^)g27re,(=r;„)C(a=„,/) ji^jg 
provides a local frequency description of component / at Xq. The instantaneous frequency should be 
interpreted in terms of the local frequency of the oscillations, |(/)j(a;o)| , and the orientation of these 
oscillations, ni{xo). The sign of (t>'i{xo) is taken so that the orientation angle is restricted from — ^ 
to |. Finally the local magnitude ai{xo) has the interpretation of local energy presence. If the image 
is actually a sum of several AM/FM/OM terms, i.e. the signal is multi-component, as is most often 
the case of many observed images, we will not be able to use this description directly, as we cannot 
separate out the L components. The problem with multi-component signals in one dimension is much 
documented [8], if not fully resolved, but generally calls for localised methods. We shall thus construct 
local monogenic descriptions of images. These descriptions extract well-behaved orientation, phase and 
amplitude functions locally, and also have excellent statistical properties. 

VI. Orientation & Monogenic Wavelets 

A. Definition 

In one dimension the analytic Morse wavelets can be constructed from the even Morse wavelets, by 
adding i times the Hilbert transform of the original even function, to the even wavelet. The even Morse 
wavelets are invariant to sign changes, or direction, whilst the odd Morse wavelets, are naturally odd 
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functions. We construct a monogenic version of isotropic wavelets in two dimensions, based on the real 
isotropic Morse wavelets, using the Riesz transform. The monogenic wavelet transform locally defines a 
phase and an orientation for a real image at each spatial and spatial scale point, similarly to the monogenic 
signal. The monogenic wavelets are most easily represented as quaternion-valued functions defined for 
each n via 

^^1+) {x) = ^(^) {x) + {x) + (x) = ^(^) (x) + (x) , (33) 

where the Fourier transform of the real part of the monogenic wavelet is given by (23) and the additional 
two real functions are defined in the Fourier domain as Riesz transforms of the isotropic wavelet function 
via 

(/) = -j^^(2vr/)'e-(2./)'"L:'- (2(2vr/)-) , . = 1, 2. (34) 

Note that we consider ipn^ (x) = i^n^ (x) + jtpn^ {x) , as a single object, and this has the same norm 
as tpn^ (x). In subsequent analysis we fix {l,m), and henceforth suppress their value, for notational 
convenience. For the case {l,m) = (8, 3) plots of the Riesz transforms in the xi direction of the isotropic 
wavelets n = 0, 1 are given in the spatial domain (see Figure 2) as well as the spatial frequency domain, 
where their modulus is plotted. Note that the real component is an isotropic wavelet (like the even 
wavelet in one dimension) and the two components are odd in the xi and X2 direction respectively. 
We define the translated, rotated and dilated wavelet 7p^^^{x), as the appropriate sum of translating, 
rotating and dilating its real valued components. The continuous wavelet transform of an image g{-) 
with respect to either the components of, or with respect to the full quaternionic wavelet, is defined 
as Wn\$;g) = f (fxg{x)ip^^* {x). The associated scalogram is S^\i;g) = \w'k\i;g)\^. The wavelet 
transform of image g{-) is then given via either the spatial domain, or spatial frequency domain, in terms 
of C, = [a, ^, /b]^, where /b is the Fourier variable for h via 

w^+^{i;g) = wl^Ht9)-iw^nHtg)-jw^nHtg), (35) 

wi+HC,9) = [l-kcos{^-e) + sm{^-e)]wi'^{C;g). (36) 

Thus the filtering carried out in the Fourier domain, can be understood, by looking at the wavelet function 
in the Fourier domain: 



^^^^nif) = [1 - cos((^ -9) + sin((^ - 9)] if) , (37) 



with modulus 



=2(l + sin(<^-0)) ^'I'Jf) 



(38) 
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2-n„n., ^^^^ 



Hence in terms of b and a the monogenic wavelet is filtering the image identically to the isotropic real 
wavelet, whilst the term sin{(j) — 6) is positioning the image, in orientation, in relation to the wavelet. This 
representation clarifies that the rotation by angle 9 of the wavelet function, as repositioning the axis of 
analysis, by a rotation of 6. Note that the modulus of the real isotropic wavelet is invariant with respect 

to e. 

Finally consider the joint structure of the N wavelets. We obtain that (see Appendix A) 

(V'i'\V'i'^) = 0, (vif.V'S) = 0, 1 = 1,2 

The multiple Morse wavelets thus form an orthogonal system, and this will have implications for their 
usage when performing estimation of local characteristics of real images. The total energy of the image 
using the nth wavelet is given by 

= Si'^it9) + Si'H$;g). (40) 
Also note that, as V'^'^H") radially symmetric, 

w^r^H^;9) = cos{9)wi!\^o;9)+MO)wi?i$o;9), 

w^^\i;g) = -sm{9)wl^\$o;g)+cos{e)wl?\$o;9), (41) 

where = (oiO, 6). Thus the wavelet transform needs only be calculated for one orientation, and can 
then be formed for any orientation by judicious recombination. We consider now the analysis of typical 
image features with the monogenic Morse wavelets. 

VII. The Monogenic Wavelet Transform of Discontinuities 

We consider both point and line discontinuities. The idealised version of a point discontinuity at xq 
corresponds to: gs^i{x) = ai{x)6{x — xq), where ai{x) is assumed to be a well-behaved function at 
point X = Xq. This singularity is characterised by xq, its location, and ai{xQ), the amplitude of the 
location. The wavelet transform of this object is w^~^\^; g) = ai(iCo)V''^^^* {f-e {^o — /O'- This is 

clearly maximum near 6 = a^o where it has a modulus square of of (tco) V'^ (0) /a?, and hence point 
singularities can be located by finding maxima in h. Furthermore g)^ has no dependence on 9, 

and as the magnitude of the wavelet at the origin is known, ai(-) can be determined. A one-dimensional 
singularity is modelled as 

9s,2{x) = a2{xi)5 {cos{62)xi + sm.{92)x2 - c) . (42) 



Febmaiy 2, 2008 



DRAFT 



STATISTICS SECTION TECHNICAL REPORT TR-05-02 



18 



The Une X2 = c {62 = vr/2) modulated by the value of a2{xi) is permitted, however we do not permit 
the line xi = c {62 = 0), as this would lead to an image of infinite energy. We will anticipate a further 
constraint on 62, that is — ^ < ^2 < X' which we will explain at the end of this section. Assume that 
02(xi) is a symmetric function around a maximum at x\ = xi^max- We characterise the structure of 
gs,2{-), using the wavelet transform. The wavelet transform using the isotropic wavelet only, noting that 

Tpn\x) = Tpn^rix), (whcrC tpnlix) Vx > Xr) is 

'Wn\^„9s,2) = X!^ a2{xi)/a'4!^nl (^ay^ (xi - + (ccsc(02) - cot(6'2)xi - 62)^) dxi. This will be 
large for values of b such that (xi_max — ^i)^ + (ccsc(02) — cot(6'2)xi^max — ^2)^ < Xr/a. Similar 
results hold for Wn '^ 93,2), as the Riesz transforms roughly have the same spatial support as the original 
wavelet. Hence we identify the location of the singularity for any fixed value a as 6 = (xi max, ccsc(02) — 
cot(02)a^i,max)i as maxima in the modulus of the wavelet transform using the monogenic wavelet. The 
orientation 62 will thus visually be apparent from the b plane but can also be characterised at a fixed 
point ^. The Fourier transform of (42) is 

r (-F\ - ^2(/l - COt(6'2)/2) -j27r(/i cos(6>2)+/2 sin(e2))c„-.727r(- sin(e2)/i+cos(e2)/2)ccot(e2) 

~ |sin(02)| 

where A2{-) is the Fourier Transform of a2(-). The wavelet transform of this image is Wn\Co] 9s,2) = 
jfy-a'^n\afi))Gs 2{f), where = [a,0,fb\. Thus we find that the Fourier transform of the rotated 
wavelet is given by wi^^(C; 5^,2) = j^^^^^^^^^^J^^^^^^o^n V/fe)Gs,2(/), with the obvious extension for 

(2) 

Wn {C J 9 3,2)- Now we choose to evaluate the wavelet transform with 62 = 0^ and define $,2 = [(^,^2,^] . 



sin(6')| y„oo fe,b 



00 

'OO 



1 r ,-4Mei2./...[(.-.5).-c]^^(^2 ^^^^^) 

\sm[9)\ 7„oo fe,b 



where the last equation defines j42(/g ^ i',$2,c), as an even function of fe^b,!, and f, = [fg^h,i, fe,b,2] = 
f-efb- Note that for fixed values of a and 9 = 62 we can find a value of b such that (r_gfe)i = c, denoted 
62- When ^ = ^2 = ^2, b2] then sin{(27r/5i ^,1) [{r^gb)i — c]} vanishes identically for all Hence 
'^n\C2'^ 93,2) = 0, whilst from equation (41) we find that the energy of the wavelet transform with ^('?)(-) 
is conserved under rotations. Thus lif'-^H^iS)! maximum at 9 = 62, and & = 62- With the correct 

(2) 

choice of rotation 9 we find that we can retrieve the angle 92 via maximising the energy of Wn {^',9s,2) 
and minimising the energy of 'Wn\^; §3^2)- If the magnitudes of the two wavelet transforms at ^ = ^0 
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are equal then 9 = Stt/A. Otherwise we take a value of 6 that maximises 

-2sm{2e)w^^\^o;gs,2)w^^H^o;9s,2). 



(43) 



This has a stationary point at 



1 

— tan 
2 



2wl^\$,o; gs,2)w^^ {^o;9s,2) 



(1)2 \ (2)2 
Wn' {^o;9s,2) -Wii' {^o;9s,2, 

which corresponds to a maximum by choosing the appropriate solution. Note that if 



(44) 



Wn\^0;9s,2) 

Wn\^o;gs,2) 



{2)^t 

wk {^o;gs,2 



we choose the solution — f < < f that corresponds to a maximum whilst if 



> 



Wn\io;9s,2) 



we take ^ < < ^, or— | < 6 < — ^. Thus at any fixed point ^ we can find the 



orientation that would result from a line-discontinuity passing through b, by utilising the above equation, 

and this characterises local orientational structure. Note also that at ^ = ^2 this evaluates the wavelet 

2 



transform at values of b corresponding to the discontinuity, 



parameter 0max,n should be combined with a plot of 
corresponding to an edge is present. 



'[$2'^9s,2, 



Wn '{^2'^9s,2, 
2 



will be large. A plot of 



to verify that local rapid variation 



< 



VIII. The monogenic Wavelet Transform of AM/FM/OM Images 

Consider analysis of images of the form given by equation (29). For more general classes of images, 
i.e. such as images that are constrained to be positive we may add a constant term to the model, but 
as noted by [14], the wavelet transform is a zero-mean filter, and so this makes no difference to the 
subsequent analysis. We then find that 

1 



ai{b) 



,2MUb)+^[{b)m{b)ix~b)) 



ai{b) cos[2^(/)/(6)]aM'(f) [a<t>'i{b)) , 



+ e 



-23iT{Ub)+'l>[{h)m{b)ix~b)) 



ai{b) 



^2jiT{Mb)+<l>[{b)Mb){x~b)) ^ ^-2j7r{Mb)+4'[{b)n,ib)(x-b)) 



ai{b) e2j'"<^'('')a^(,™> {a(t>'i{b)ni{b)) + e-^^^'''''^''^ a^'^ 



h\{b)ni{b)) 



m 



1,2. 



Hence it follows from equation (41) that w^+\i;ci) = ai{b)a'i!^nHa(t>'i{b)) e^'^''"M^\ with e^^ = 
i cos {r]i{b) — 6) +j sin {rii{b) — 6). This is the localised analogue of equation (31). For multi-component 
images, we may be able to separate the dominant component out, similarly to ridge analysis based on 
complex wavelets [14]; this requiring the assumption ai{b)a'^n'' {acpUb)) >> am{b)a'^n^ (acp'Ab)) V / 7^ 
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m, at all ^ considered. Furthermore, the modulus of the wavelet transform is 

w(-+\$;ci)f = a2(6)a2*(f)2 (a<^'^(b)) , (45) 

and hence the wavelet transform of q(-) is locally maximal on the curve given by 
n{a,9,b) = [{a,b): a^^^) = /mix} , 

where /]^2x is given by (24). This defines the monogenic 
wavelet ridges [14] of an AM/FM/OM images. At any point on this ridge, the local orientation may 
be computed. Ridge analysis is based on the fact that not all information of the redundant wavelet 
transform representation needs to be considered to characterise the image: only the ridge itself. As the 
ridge definition does not depend on the angle 9, we need not carry out the transform for all these values; 
a computational advantage to using directional wavelets. On the ridge we characterise the oscillatory 
components locally as 

ui{b) = tan-^ 1 =m{b)-e, (46) 



w 



^'^'^ = [ ^m^) j' ^ ^ 

We have constrained — f < i^; < f and — ^ < (/) < ^ by the choice of sign for the tan~^(-). 

IX. Digital Implementation 

To preserve the exact monogenic structure we implement the wavelet transform from the Fourier 
domain, calculating the IDFT, thus making the algorithm of order 0{NiN2 log(A^i) log(A^2))- We consider 
the maximum and minimum scales that can be resolved - the range of the angle 6 G (0, 2tt) and 
b G (0,A'^iAi) X (0,A''2A2). As the real two-dimensional even wavelet "^n\-) is built from a real 

(n) (n) 

one-dimensional wavelet corresponding to a band-pass filter, there exist frequencies /{ and such 
that 

*^f^(/)^OV/: l/l^ (/f ,/2^"^). (49) 

Note that the DFT of observed image g{-, •) is periodic by construction, and that the standard assumption 
corresponds to /2) = for all frequencies not in the Nyquist band. Thus to perform the implemen- 

tation we consider only scales a such that ^'^f^ (a\//i + /I ) ~ V / : /i > l/(2Ai), /2 > 1/(2A2). 



This necessitates a > am\„ = 2fi"^— ^i^2=. As a increases in magnitude, the wavelet becomes more 
peaked in the frequency domain, and to ensure the wavelet covers at least M frequency points we 
constrain a < a^ax = jj min (iViAi, iV2A2) /a"'' - f^^ 
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X. Statistical Properties 

Consider estimation of features present in an image immersed in white noise where the image is 
collected in a regular grid consisting of xi = siAi, si = 0, ...,A'^i — 1 and X2 = 82^2, S2 = 
0, . . . , N2 — 1. We model the observed image y{xi,X2) as ?/(xi, X2) = g{xi,X2)+£si,s2 - The noise {£51,82} 
is modelled as isotropically Gaussian and white. It is assumed that E (esi,^^ = ^i^d ^ (£si,s2) ^ui,u2) = 
'^si.ui'5s2,«2<^e • The wavelet transform of the noise will also be Gaussian, as it corresponds to a sum of 
jointly normal variables. To give the distribution of the wavelet transform, we calculate its first and second 
order structure at a fixed ^. The wavelet transform is a linear operation and Wn\^;y) = Wn\^;g) + 
Wn\^;e). It follows that E (^Wn\$,; y)j = Wn\$,;g), and we can determine the second order structure 
of the estimators from the distribution of the noise. 

We recast the full wavelet transform of the noise as a vector with real valued entries, i(;„(^;e) = 

n = 0, . . . , N — 1. In Appendix B, we find with the additional 



assumption of 



that 



mm 



^ ^ ' > max ^ (50) 



2Ai 2A2y n=0,...,N-l 



'^10 

Wn{^;e)=Af3{03,a^,V), V= i 



(51) 



We perform estimation using the multiple orthogonal wavelets. Any estimator of local signal properties 
needs to be smoothed, or averaged to obtain a low variance. The wavelet transform using any of the 
specified wavelet functions averages the data across a window in space and spatial frequency, where the 
width of the region depends on the wavelet chosen, and in our case is characterised by the radial Morse 
region V, and the parameters (/3 = / + 7 = m). Thomson [1] suggested forming estimates of local 
properties by averaging local energy estimates using several orthogonal wavelets/functions. This usage 
explicitly reduces the variability of the estimates with a clearly specified averaging region - V. Coherent 
behaviour over V is re-enforced across wavelet estimates, but the noisy uncorrelated behaviour should 
cancel. The bias inherent in the averaging is characterised by the eigenvalues square of the localisation 
operator. In Appendix B we show that E (^Wnl\^', e)wn2\^,'i = ^'^'^'5ni,n2) ^nd thus it;„i(^;e) 
is uncorrelated with i(;„2(^;e) unless rii = n2, that combined with the assumption of Gaussian errors 
corresponds to independence. We define averages of the wavelet transform and the scalogram that will 
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be used as a basis for calculating estimators of other quantities as 



n=0 n=0 

with / = e, 1,2,+, and finally as a measure of covariation we define for h, I2 = e, 1,2,+, and 
images g2{-), c'-^^'^^ht gii'), 92i-)) = ^ En=o^ ^^"'^(^i 5i)'f^n'''*(^; 52)- We define the estimators 

u)(0(^;c,) = iu(0(^;y) and S^^\^;g) = S^^\ty), for / = e,l,2,+, as well as w(^^^^) i^; gi) = 
^('i,'2)^^_ g^(^.^^ g-^(^.^'^_ jhe Gaussian assumptions on e then give e) = e),w^^\^; e),w^'^\^,; e)] 

M (O3, CT^V) , where = al/N. For most quantities we would intuitively expect to see a reduction of 
1/A^ in their variances. When estimating the energy of the image at point ^ we consider 



d 



N 

n=0 



2 



^^"-^ ; 5) + 5^+) ; e) + _ ^ ) ; e) (54) 

(55) 



n=0 

Up to order ct^, we find, with the additional assumption of the localised behaviour of g{) coherent across 
the n wavelets, 



2 

Hence the variance of the energy estimate decreases O (■^) . 



S^"\^;g) + Us^'\^'9)+S^'\^;g) 



(56) 



A. Distribution of Estimators 

We estimate the orientation of the line discontinuity in section VII by maximising the difference 
between the energy of the second and first components. Each wavelet indexed by n satisfies equation 
(44) and thus writing the equations in terms of tan(0max,n) = 02, we may sum over the equations to 
find that 

n K -if 2C^''^\t9s,2,gs,2)) \ 

(72 = — tan 

We form estimate 



2 W'\io;gs,2)-S^'\^o;gs,2) 



1 -1 f 2w(i^)w(i^){^;gs,2) 



2 \S^'K^o;9s,2)-S(^){^o;gs,2)J 
Let 'Wn\$,;e) = cr^Wn\, which entails that w^^\^;e) = a^wT* for I = e, 1,2,+, and expand the above 
expression 0max(s's,2) = 6*2 + (ye562 + 0{a1). Note that E {602) = making the estimator up to order 
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o"^ unbiased, and the estimator has variance 



Var 



^1 S^'\gs,2) + S^'\gs,2) 

N 2 



-(1), 



+ 4 



7=7(1.2)/ . 
C {gs,2,9s,2] 



+ 0{a! 



(57) 



S {ys,2 

Thus, using multiple wavelets leads to a variance reduction. 

For AM/FM/OM signals we define the estimator for the orientation angle of the unit quaternion as 



n{i]ci) = i^iitVi) = tan ^ 



that up to order is 



(2) 



(58) 



+ 0[al) 



(59) 



as w^^\$,Q]Ci) = ai(6) sin(27r(?l);(&)) cos(f/)-^ ^ (a I V<^i(6)|) and a similar results holds for the 



second component we have tan(i/i) 



'AT 



zero the estimator is unbiased and it has variance 



. As the wavelet transform of the noise has expectation 



Var[I?z(^;Q)] 



(2)2/ 



2N 



'S^'\i;ci) + S'^'{i;ci) 
al (E^n.(a|V,^,(b)l))' . 



o(2), 



0". 

2Af 



(60) 



(61) 



2iv2 E^i(«|v<^Kb)l) 

Using multiple wavelets leads to variance reduction. To estimate the phase we only use a single wavelet, 
the n = 0. Due to the orthogonality relations, the wavelet filters in the Fourier domain cannot be strictly 
positive for all frequencies, and thus forn > there are induced variations in the phase estimate whenever 
the wavelet filter changes sign. The Morse wavelet estimate of the phase still profits from the wavelet's 
good radial localisation. 



(62) 



27r 



w 



e,0 



,(i),„(i) 



(2),„(2), 



As the expected value of the wavelet transform of noise is zero, the estimator is thus unbiased, and the 
variance of the phase estimator is 



Var 



al 5i+)(C;Q)-i5^^^(g;Q 

:(+) 



1 o{e)/ 



(27r)25j+)(^;Q) 



(27r)2si+)(^;q) 



1 



■cos2(27r(/)/(6)) 



(63) 



When considering larger scales, the wavelets are averaging across a lot of sample points, and the variance 



of the phase estimate decreases. The amplitude is estimated as af{b) 



a2*(=)2(a(/);(b)) ' 
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Fig. 3. The local energy of signal 1 on a dB scale using three wavelets (far left), or one wavelet (second from left), a = 1.4. 
The deviation of the estimated orientation from the true orientation at scale a = 1.66 using dB scale, using three wavelets 
(second from right) and one wavelet (far right). For the region where the signal has presence at those spatial and scale points, 
the estimate using three wavelets is less noisy. 



XI. Examples 

Consider a collection of singularities observed in noise: gi{x) = '^g2j{x) + criex where gii{x) = 

, i°ijv, + H"ll ' = i|^_r45jv, + i^i^jv, + ir"ir ^^"^^""^ = '"^^ cos(7r/3) - X2 sin(7r/3) - + \ 

" U 4 ^^2j II II [64 2 64 ^^2j I 

gi4:{x) = |a;i cos(7r/9) — X2 sin(7r/9) — ||-/Vi + W , and we take ai = 0.2. The first two singularities 
of this signal are point singularities and the latter two are line singularities. See Figure 3 for a plot of the 
scalogram of the observed image at scale a = 1.4, corresponding to radial frequencies of 0.17. It can be 
clearly made out that the averaged estimate of the local energy is a great deal more robust to the noise. 
Signal 2 is a multi component AM/FM/OM signal given by g2{x) = g2i{x) + 522(2;) + a^ex, where 
g2i[x) = 1.21 (xi < iVi/2) cos (itt x 0.087(2!^ + ts)) , U = xi cos(7?0 + X2 sin(r?,), / = 1, 2, r/i = 
-f + ^^tIiVT^' 922{x) = 0.8 cos (o.OSvr + ^2)) , ??2 = f , and Ni = N2 = 128. We consider 

estimating its orientation at a scale where the more rapid sinusoid is present near the left-hand side at 
those frequencies, and find that our orientation estimate is substantively less noise when using multiple 
wavelets, as is confirmed by Figure 3. Clearly using the multiple wavelets is substantively decreasing the 
variability of the estimator. 

XII. Conclusions 

We show that the multiple monogenic Morse wavelets hold great potential for digital image processing 
and analysis. The monogenic Morse wavelets are the natural two-dimensional extension of the analytic 
Morse wavelets, and are the eigenfunctions of a two-dimensional, nonseparable, localisation operator. 
They form an orthogonal system, where the orthogonality establishes the statistical properties of Gaussian 
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noise. By averaging across wavelets, estimators of local properties of the signal achieve reduced variability. 
The monogenic properties of the wavelets form a natural framework for determining local phase and 
orientation properties. This framework explicitly parameterises the local orientation of any variational 
structure via a unit quaternion, and the localised analysis considers radial structures, thus using the 
natural metric of Cartesian distances in the spatial domain. 
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APPENDIX 

A: Orthogonality Relations 
We know that if we choose c = {21 + l)/m — 1 [3, p. 2663] then 




^-2(2.1/1)'" (2(27r I/I)-) L^, (2(27r |/|)™) df 




ni ,n2 ■ 



(64) 



Consider the two-dimensional integral of two dimensional wavelets 




IT 



1 




An,;l,m{27Tf) 



'e-(2-/)'"L^; (2(2^/)™) A„,,,,„,(2^/)'e-(2-/)'"L5^; (2(27r/)-) dh dh 



poo 

2 / A„,,,,^A„,,;,^(2vr/)2'+ie-2(2./)-^c'^ (2(2vr/)-) < (2(27r/)™) df 
Jo 




'711,712- 



(65) 
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Also note that 

(•OO /"OO 



(^^V'S) = - / / A.,;^™(27r/)'e-(2-/)'"<(2(27r/)^ 



OO ^ — OO 

'2 



A,;i,^(27r/)'e-(2-/)"L^^^ (2(27r/)-) d/i d/2, / = 1,2. 



/f + /I 

This impUes that (t/;^^,^ , V^is^ ) + (V'n? , -^n} ) = {ipnr , il^n2 > = ^rn ,n2 , and as ^-^f ^ (/) is radially symmetric 
we may deduce that 

(V^«,V'li)) = (^i^,),Vl'.)>, (66) 

and thus 

{¥n'^,4!^) = lSn„n2. (67) 

Finally, note that 

^ J —OO J —OO 

A„,;;,^(27r/)'e-(2./)™^c'^ (2(27r/)™) (-i)^^^?! 

J I + /2 

= 0, (68) 
due to the integral of an odd function over a symmetric region being zero. Similarly 

(V'if,^i^)) = 0, (^W,V^(2))=0. (69) 

B: Calculation of Statistical Properties 
Define the discrete Fourier transform of the noise ex^,x2^ 

£{fij2) = E E ex„x.e-2i-(^^^+^^^^). (70) 

Xi=0 X2=0 

As the wavelet transform at any angle 9 can be formed from linear combinations of the wavelet transform 
at = 0, in the way outlined in section VI, we need only calculate the properties at = 0. We have: 



2 2 ^i^'^ 

-Ni l3=-N[ U=-Ni 



" V V V , S, , p'^Mblill-l3) + b2il2-U)) 

]vr2 Tvr2 a2 /\2 '^liM^hM^ 



1 1 



2Ai 2A2 



aV?ai / / ( a.//2 + /f ) ( aJf^ + /| ) d/i df2 (71) 



0"^5ni,n2, (72) 
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where the last Une follows from Eq. (65). Similarly 

Covf<)(|;e),w;W(^;e)' 



2Ai 2A9 



— 2 <:"rii,n2- 

following from Eq. (67), and likewise 

Cov(u;i2)(^;e),^(2)(^;e)) = ^a^J, 

Furthermore, 

Cov(^.(f(^;6),^.W(^;e) 
= E (w^^X<^^ 0' 61 Ai, 62A2; e)w^^>{a, 0, 61 Ai, 62A2; e 



2 - e ^ni,n2 • 



2Ai 2Ao 



0, 



^ 



following from Eq. (68), likewise Gov [^fWrii e), u^na (^; e) ) = according to (69), and finally 

Cov(<)(^;6),t«(2)(^;e)^ 
= E 0, 61 Ai, 62A2; e)w^;^> (a, 0, 61 Ai, 62 A2; e)^ 

= 0, 

following from Eq. (69). This completes the covariance calculations for the distribution. 
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